function [optima, path, budget] = LocalSearch_BoxConstraint(obj,bounds,starts,step,step_threshold)
%LocalSearch_BoxConstraint: the local search algorithm for box constrainted
%   continueous minimization problems
%neighboring points are defined as the points where one of the dimension moves a step
% INPUTS:
%   obj: the objective function handle, f = obj(x)
%   bounds: 2 -by- d matrix, [lb;ub]
%   starts: n -by- (1+d) matrix, the start points, [objective value, x1, x2 ... xd]
%   step: the value of one dimension of the neighbor points changes
%   step_threshold: stop when step<=step_threshold and the local optima is
%       the current point
% OUTPUTS:
%   optima: n -by- (1+d) matrix, the final optimal solutions, [objective value, x1, x2 ... xd]
%   path: n -by- 1 cells, the path of the current best point
%   budget: n -by- 1 vector, the total budget used from each start
%       points (exclude start points)
[n,d] = size(starts); % the number of start points
d = d-1; % the dimension of the decision vector
optima = zeros(n,d+1);
path = cell(n,1);
budget = zeros(n,1);
for i = 1:n
    currentbest = starts(i,:);
    path{i} = currentbest;
    while 1
        currentbest_changed = 0;
        for j = 1:d
            % generate neighboring points
            neighbors = [currentbest(2:end);currentbest(2:end)];
            neighbors(1,j) = neighbors(1,j)-step;
            neighbors(2,j) = neighbors(2,j)+step;
            neighbors(neighbors(:,j)>bounds(2,j),:) = [];
            neighbors(neighbors(:,j)<bounds(1,j),:) = [];
            NumNeighbor = size(neighbors,1);
            % objective values
            fvalue = obj(neighbors);
            budget(i) = budget(i)+NumNeighbor;
            % find the best neighboring point
            [~,best_temp] = min([currentbest(1);fvalue]);
            if best_temp==1
                path{i} = [path{i};currentbest];
            else
                currentbest_changed = 1;
                currentbest = [fvalue(best_temp-1),neighbors(best_temp-1,:)];
                path{i} = [path{i};currentbest];
            end
        end
        if currentbest_changed==0
            if step<=step_threshold
                break
            end
            step = step/2;
        end
    end
    optima(i,:) = currentbest;
end
% remove repeated optima
i = 1;
while i<=size(optima,1)
    dis = sqrt(sum((repmat(optima(i,2:end),[i-1,1]) - optima(1:(i-1),2:end)).^2,2));
    if min(dis)<step_threshold
        optima(i,:) = [];
    else
        i = i+1;
    end
end

end

